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1  INTRODUCTION 

Seaskimming  missiles  fly  as  long  as  possible  below  the  radar  horizon  of 
the  target  ship  to  reduce  the  detectability.  Therefore  the  missile  has 
to  fly  as  low  as  possible  at  constant  height  above  the  mean  sea  level 
without  ditching  the  water. 

To  fly  at  a  constant  altitude  the  missile  has  to  be  equipped  with  an 
altimeter  to  measure  the  missile  height  above  the  sea.  The  altitude 
control  law  imbedded  in  the  guidance  computer  controls  the  actual 
distance  between  the  missile  and  the  water  surface.  In  figure  1  the 
situation  above  the  waves  is  depicted. 

To  simulate  the  flight  above  the  sea  waves,  these  waves  have  to  be 
modelled  and  simulated.  The  sea  wave  motion  can  be  described  by  a 
stochastic  process  and  is  determined  by  i.a.  wind  speed  and  depth  of  the 
water . 

In  chapter  2  a  description  is  given  of  the  power  spectral  density  (PSD/ 
function,  which  describes  the  power  distribution  of  the  sea  waves  as 
function  of  frequency  for  a  stationary  and  moving  observer.  In  chapter  3 
some  methods  are  described  for  the  generation  of  sea  wave  elevations.  In 
chapter  4  the  influence  of  sea  wave  motion  on  missile  performance  is 
analysed.  In  chapter  5  some  conclusions  are  given. 

This  study  has  been  carried  out  under  assignment  nr.  A84KM080  "Missile 
simulations".  Results  have  been  used  in  various  studies  for  instance  the 
simulations  of  the  threat  for  the  NATO  Anti  Air  Warfare  System  (NAAWS). 
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2  POWER  SPECTRAL  DENSITY  FUNCTION  FOR  STATIONARY  AND  MOVING 

OBSERVER 

The  sea  wave  motion  can  be  described  by  a  stochastic  process.  This 
process  is  specified  by  a  power  spectral  density  (PSD)  function  and  a 
probability  density  function  (PDF). 

We  will  now  consider  the  sea  wave  PSD  function  with  normal  distribution. 
The  state  of  the  sea  is  described  by  the  sea  state  number,  see  table  1. 
The  sea  state  is  defined  by  a  parameter  known  as  the  significant  wave 
height  which  is  the  average  of  the  highest  third  of  all  waves 

recorded  in  a  specified  interval,  see  ref.  1. 

The  PSD  of  the  sea  wave  motion  is  given  by  the  Pierson-Moskowitz  (PM) 
spectrum  and  is  a  standard  accepted  description  of  this  stochastic 
process.  It  is  a  one  dimensional  spectrum  for  a  stationary  (non-moving) 
observer  for  the  case  of  a  fully  developed  sea  in  deep  water  and  is 
given  by  : 


-3.36 

0.84  H3rd2-u.4 

$  (w)  - - e 

wv  '  S 


The  standard  deviation  a  equals  : 

^3rd  -00 

a - J  (  /  4>  -dw  ) 

4  0 


(1) 


(2) 


The  significant  wave  height  for  sea  state  (SS)  2  is  0.6  meter  and  for 
SS  5  equals  4.0  meter  according  to  table  1.  For  simulation  these  two 
values  are  most  commonly  used.  The  PSD  functions  for  these  two  states 
are  given  in  figure  2  and  3. 

The  above  given  spectrum  (equation  1)  is  for  a  non-moving  observer.  To 
account  for  missile  velocity  the  spectrum  has  to  be  transformed  to  a 
spectrum  for  a  moving  observer,  because  the  altimeter  moves  with  the 
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missile  at  a  speed  V  as  it  senses  the  height  above  the  sea  waves.  The 
frequency  for  a  moving  observer  is  shifted  according  to  : 


GJ  -  (D  +  W 

m 


w-  (1  +  -  ) 

C„ 


with  Cw  the  velocity  of  the  wave  propagation  and  U  is  the  horizontal 
component  of  the  missile  speed,  towards  the  direction  of  the  wave 
propagation.  The  velocity  of  wave  propagation  is  : 


with  g  the  gravity  acceleration.  So  the  frequency  for  a  moving  observer 


u  -  w  (1  + - w) 

m 


The  spectrum  for  a  moving  observer  then  becomes 


$  ( )  — 


2U 

1  + - w 

g 


The  derivation  of  this  spectrum  is  given  in  appendix  A. 


To  generate  sea  waves  the  process  as  indicated  in  figure  4  is  used. 

White  normal  noise  is  filtered  to  get  the  sea  wave  elevations.  The  power 
spectral  density  function  of  this  process  is  given  by  : 

*w  (»)  -  |F(jw)|2-*N  (7) 


is  the  power  spectral  density  function  of  white  noise. 
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The  problem  is  now  to  find  a  filter  that  matches  with  the  PM  spectrum. 
In  ref.  2  a  frequency  response  F(jw)  is  found  according  to  : 

b-(jw)4 

F(ju)  -  - —  (8) 

[(jw)2  +  J 2-a-jw  +  a2]3 


The  parameters  a  and  b  are  defined  as  : 


a 


max 


4 

/  J  2 


b  = 


a2  9 

_  j  ^  “max 

0.385  SN 


(9) 


where  $  is  the  maximum  value  of  the  PSD  function  of  the  sea  waves 
wmax 

for  a  non-moving  observer  and  is  the  frequency  at  which  the  PSD 

function  for  a  non-moving  observer  reaches  its  maximum  value.  In 
appendix  B  the  derivation  of  is  given. 

A  filter  which  gives  a  better  match  of  the  theoretical  PSD  (PM- spectrum) 
has  a  frequency  response  F(jw)  according  to  : 

b-(jw)4 

F(j»)  -  - —  (10) 

(jw  +  a)b 


where  the  parameters  a  and  b  are  now  defined  as 


a  -  cj 


®max  /  J  ^ 


b  -  1.53V  •  J  (  — ^ax  ) 

®max  * 


9 


“max 


(11) 


The  filter  as  given  in  equation  (10)  will  be  named  FEL  filter. 

In  the  figures  5  and  6  the  results  of  spectra  calculations  are  given  for 
SS  5  with  an  observer  speed  of  340  and  680  m/s,  respectively.  In  these 
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figures  the  spectra  according  to  Pierson-Moskowitz  (PM),  to  ref.  2  and 
to  the  FEL  filter  are  given.  The  PM-spectrum  is  calculated  with 
equation  (1).  The  other  spectra  are  calculated  with  equation  (7).  The 
spectrum  made  with  the  FEL  filter  is  closer  to  the  PM  spectrum. 

For  sea  state  5  the  variance  of  the  signal  is  one.  The  variance  is  equal 
to  the  integral  of  the  PSD  function.  The  filter  from  ref.  2  gives  a 

O 

variance  for  SS  5  equals  to  0.66  m  and  the  FEL  filter  gives  a  variance 
value  of  0.93  m  .  A  criteria  to  determine  how  good  different  solutions 
match  is  the  weighted  integral  squared  error  (ISE)  criterion.  The 
equation  for  this  criterion  reads  : 

GO 

ISE  -  f  WF(u>)  ■  l*u(u>)  -  fcpnCuO]2-^  (12) 

0 

WF(w)  is  the  weight  function. 

The  ISE  is  dependent  of  observer  speed.  The  ISE  at  340  m/s  for  ref.  2 
is  approximately  5  times  as  large  as  the  ISE  for  the  FEL  filter  if  the 
weighting  function  is  one.  If  the  weighting  function  equals  six  from 
zero  to  w$max  and  one  for  the  other  part  of  the  frequency  spectrum  than 
the  ISE  for  both  filters  is  the  same. 

The  FEL  filter  gives  a  better  match  according  to  the  integral  square 
error  criterion  and  will  therefore  be  used  for  the  generation  of  the  sea 
wave  elevations. 
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3  METHODS  FOR  GENERATION  OF  SEA  WAVE  ELEVATIONS 

Sea  wave  elevations  or  sea  wave  heights  are  calculated  with  the  process 
described  in  figure  4,  but  there  are  different  methods  for  doing  this. 
Two  methods  are  described  here,  we  call  the  first  method  generating  sea 
wave  elevations  on-line  and  the  second  method  generating  sea  wave 
elevations  off-line.  In  the  off-line  method  an  on-line  correction  for 
missile  speed  is  applied. 

The  first  method,  which  was  formerly  used,  is  based  on  the  calculation 
of  a  batch  of  sea  wave  elevations  as  seen  by  a  moving  observer  for  one 
second  with  a  sample  time  of  0.01  sec,  so  this  batch  contains  100 
samples.  This  is  done  every  second  again  with  newly  calculated  filter 
coefficients  because  of  possible  fluctuating  missile  speed.  This  method 
has  several  drawbacks.  Firstly  it  is  very  time  consuming  because  every 
second  again  the  filter  has  to  be  settled  before  the  filter  can  generate 
sea  wave  elevations.  Due  to  the  large  time  constants  involved  1100 
samples  have  to  be  generated  of  which  only  100  are  used.  Secondly  it  is 
less  accurate.  The  PSD  function  for  this  method,  calculated  from  an 
autocorrelation  function,  is  given  in  figure  7  (observer  speed  is 
340  m/s).  The  data  file  contained  20000  samples.  As  can  be  seen  low 
frequencies  are  introduced.  In  appendix  C  the  derivation  of  the  PSD 
function  from  N  data  points  is  described. 

The  second  method  is  based  on  the  generation  of  sea  wave  elevations  for 
one  speed  only,  with  a  correction  afterwards  for  the  influence  of  the 
actual  speed,  if  the  actual  speed  and  the  speed  used  in  the  generation 
differ . 

In  figure  8  the  spectrum  is  given,  also  computed  from  an  autocorrelation 
function  (observer  speed  is  340  m/s).  If  figure  7  and  8  are  compared,  it 
is  clear  that  method  2  is  more  accurate,  because  less  low  frequencies 
are  introduced. 

The  term  "travelled  distance"  is  introduced  to  correct  for  missile 
speed.  This  correction  is  necessary  if  this  speed  differs  from  the  one 
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used  for  the  generation  off-line.  The  travelled  distance  is  a  function 
of  missile  speed  and  actual  time  of  flight. 

To  illustrate  this  method  an  example  is  given.  With  the  off-line  method 
sea  wave  elevations  are  generated  for  sea  state  2  and  a  missile  speed  of 
170  m/s.  Thu  sampling  time  is  0.01  seconds.  To  prevent  aliasing  the 
sampling  frequency  must  be  larger  then  2  times  the  highest  frequency  in 
the  spectrum.  Every  sample  corresponds  with  a  distance  of  1.7  meter. 

This  means  that  the  sea  wave  elevations  are  now  only  a  function  of 
travelled  distance. 

Now  suppose  that  the  missile  is  flying  at  340  m/s.  From  the  PM  spectrum 
a  highest  frquency  of  450  rad/s  or  72  Hz  follows.  The  sample  time  must 
therefore  be  smaller  than  0.007  sec,  to  prevent  aliasing.  The  sampled 
distance  is  1.7  meter.  This  means  that  a  sample  must  be  taken  from  the 
data  file  if  the  missile  has  travelled  1.7  meter.  So  if  the  simulation 
sample  time  is  0.005  sec  then  every  simulation  step  a  sample  must  be 
taken  from  the  data  file. 

This  method  isn't  exactly  correct,  because  the  assumption  is  made  that 
the  PSD  function  is  linearly  dependent  of  missile  speed  only.  This 
simplification  means  it  is  assumed  that  the  sea  waves  which  are 
overflown  have  a  constant  (i.e.  not  time  dependent)  shape  like  terrain 
has  whereas  sea  waves  move  with  a  frequency  dependent  speed.  According 
to  these  assumption  the  transformed  frequency  is,  if  the  speed  is  x 
times  the  speed  V,  that  was  used  to  make  the  data  file  : 

x  ■  V  ■  w2 

wm  “  x'w  +  -  0-2) 

g 

instead  of 

x-V-w2 

wm  -  w  +  -  (13) 

g 

For  an  observer  speed  of  340  m/s  and  SS  5  an  error  is  made  of  0.5%  in 
shifted  frequency  if  the  missile  flies  2  times  as  fast. 
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The  power  contained  in  the  spectrum  is  the  same  in  both  cases,  only  the 
spectrum  is  not  shifted  correctly,  however  with  small  error.  To  prove 
that  the  power  does  not  change  the  next  derivation  is  given. 

If  the  PSD  for  an  observer  speed  of  V  in  the  frequency  band  [wl,w2] 
equals  $(w)  then  the  PSD  for  an  observer  speed  of  xV  in  the  frequency 
band  [xcj1,xu21  equals  $x(wx)  -  l/x'4(w/x). 

The  power  contained  in  this  frequency  band  is  : 

f  4>(w)-dw  -  f  l/x-*(u>)  •  (dxw)  (14) 

with  wx=xw  follows  : 

w=w2  wx”“x^  ux”“x^ 

f  ®<w)-dw  -  f  l/x-*(«x/x)-(dwx)  -  f  V«x)  -du>x  (15) 

w-wl  »x-wxl  ^x^x1 

This  method  was  tested  for  sea  state  2  and  5. 

The  test  consists  of  4  steps.  These  steps  are  : 

1.  Generate  a  data  file  containing  sea  wave  elevations  for  a 
missile  speed  of  340  m/s  with  a  sample  time  of  0.01  sec. 

2.  Generate  a  data  file  containing  sea  wave  elevations  for  a 
missile  speed  of  680  m/s  with  a  sample  time  of  0.01  sec. 

3.  Compute  the  PSD  function  for  the  data  file  generated  at 
step  2. 

4.  Compute  the  PSD  function  for  the  data  file  generated  at 
step  1,  but  each  step  a  sample  skipped. 

In  step  1  a  datafile  is  generated  with  a  sampling  distance  of 
3.4  meter.  In  step  2  a  datafile  is  generated  with  a  sampling  distance  of 
6.8  meter.  For  sea  state  5  a  sampling  distance  of  13  meter  is  still 
enough,  so  it  is  oversampled.  For  the  generation  of  the  PSD  in  step  4 
samples  can  therefore  be  skipped  instead  of  taking  every  sample  as 
function  of  travelled  distance. 
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Figure  9  presents  the  computed  PSD  functions  for  a  speed  of  680  m/s. 
From  this  figure  it  can  be  concluded  that  the  described  method  is 
sufficiently  accurate  for  use  in  simulation  studies. 
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4  THE  INFLUENCE  OF  SEA  WAVE  MOTION  ON  MISSILE  PERFORMANCE 

With  the  six  degree  of  freedom  (DOF)  Weapon  Analysis  and  Simulation 
Program  (WASP)  a  seaskimming  missile  is  simulated.  For  a  description  of 
WASP  see  ref.  3. 

The  missile  is  flying  at  nearly  constant  speed  at  5  meter  above  mean  sea 
level.  A  sea  state  of  5  is  assumed,  which  corresponds  to  a  significant 
wave  height  of  4  meter  or  a  standard  deviation  of  1  meter. 

Some  airframe  parameters  are  given  in  table  2.  The  time  constant  of  the 
innerloop  (acceleration  autopilot)  is  0.23  seconds  and  nearly  constant 
during  flight.  The  structure  of  the  altitude  autopilot  is  given  in 
figure  10.  The  time  constant  of  the  altitude  loop  is  in  the  order  of  2 
to  3  seconds  for  the  parameters  given  below. 

Due  to  altimeter  and  seeker  restrictions,  the  roll  angle  has  to  have  a 
constant  value  during  flight.  Therefore  the  roll  control  system  is  much 
faster  than  the  pitch  and  yaw  channels. 

Four  tests  are  performed.  For  each  test  the  autopilot  has  the  same 
innerloop  branches,  but  differs  in  feedback  filter  time  constants  for 
the  measured  height  and  filtered  height  rate. 

The  altimeter  is  assumed  to  have  no  biases  or  measurement  errors. 

The  cut-off  frequencies  of  the  filters  are  important,  because  the  high 
frequencies  in  the  spectrum  can  be  filtered  out.  For  a  missile  speed  of 
340  m/s  and  sea  state  5  the  peak  in  the  PSD  function  is  at  13.5  rad/s. 
The  time  constant  of  the  height  rate  filter  was  0.2  sec  in  batch  1  and  2 
and  0.143  sec  in  batch  3  and  4.  The  difference  between  batch  1  and  2  is 
the  filter  used  in  the  altitude  feedback  loop.  In  batch  1  the  measured 
altitude  is  not  filtered  and  in  batch  2  the  measured  altitude  is 
filtered  by  a  first  order  network  with  a  time  constant  of  0.2  seconds. 
The  same  applies  to  batch  3  and  4  but  now  with  a  filter  time  constant  of 
0.143  seconds. 

For  each  test  an  undisturbed  run  and  three  disturbed  runs  were  made.  The 
mean,  variance  and  standard  deviation  of  missile  altitude  and  angle  of 
pitch  are  presented  in  table  3.  These  values  are  calculated  from  the 
disturbed  runs  using  the  nondisturbed  runs  as  reference. 


TNO  report 


Page 

14 


To  illustrate  the  effects  of  the  choice  of  altitude  loop  parameters  on 
the  missile  performance  some  comparative  figures  are  given. 

In  figures  11  and  12  horizontal  flight  profiles  are  given  for  a  run  from 
batch  2  and  batch  3,  respectively  together  with  the  undisturbed  runs. 

Due  to  the  reduction  in  missile  mass  during  the  flight  the  seaskimming 
altitude  increases  a  little  bit  with  time. 

In  figure  13  the  input  (calculated  sea  wave  elevations)  of  the  altimeter 
is  given.  In  figure  14  the  output  of  the  altimeter  for  a  run  of  batch  2 
is  depicted.  In  figure  15  the  angle  of  pitch  for  a  run  of  batch  2  is 
given  and  in  figure  16  for  a  run  of  batch  3.  The  dashed  line  in  both 
figures  are  the  undisturbed  reference  runs. 

The  choice  of  the  altitude  loop  parameters  is  sometimes  determined  by 
contradicting  requirements.  If  the  altitude  loop  has  to  be  fast  a  lot 
of  noise  will  be  introduced  in  the  system  and  the  influence  of  the  sea 
waves  on  altitude  performance  will  be  large.  This  can  have  influence  on 
the  commanded  seaskimming  altitude. 
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5  CONCLUSIONS 

The  off-line  generation  of  a  data  file  containing  the  sea  wave 
elevations  as  a  function  of  travelled  distance  can  be  used  to  simulate 
the  sea  surface  as  sensed  by  an  altimeter  aboard  a  missile. 

The  sampling  time  must  be  correctly  chosen.  This  is  dependent  of  missile 
speed  and  sampling  distance  in  the  off-line  generated  data  file. 

The  choice  of  altimeter  loop  parameters  can  have  much  influence  on  the 
missile  vertical  flight  profile.  This  choice  is  sometimes  determined  by 
contradicting  requirements. 
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batch 


mean 

variance 

std  deviation 

(meter) 

(meter2 ) 

(meter) 

0.003 

0.00081 

0.028 

0.001 

0.00026 

0.016 

0.005 

0.00121 

0.035 

0.002 

0.00050 

0.022 

Table  3a  :  Some  results  of  the  influence  of  sea  wave  motion  on  the 
missile  altitude. 


batch 


mean 

(deg) 

variance 

(deg2) 

std  deviation 

(deg) 

-0.018 

0.408 

0.63 

-0.003 

0.079 

0.28 

-0.028 

0.647 

0.80 

-0.008 

0.195 

0.44 

Table  3b:  Some  results  of  the  influence  of  sea  wave  motion  on  the 

missile  angle  of  pitch. 


Figure  4  :  Proces  for  the  generation  of 

sea  wave  elevations 
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Fig.  6  :  Some  PSD  functions  of  sea  wav 
for  SS  5  (moving  observer  V  =  680  m/s 
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Fig.  13  :  Input  signal  of  altimeter 


Range  (km) 


Measured  height  above  the  sea  waves  (m) 
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PSD  FOR  A  MOVING  OBSERVER 

The  power  contained  in  the  frequency  band  [0  ,  cl>]  of  the  stationary 
observer  must  be  the  same  as  the  power  contained  in  the  frequency  band 
[0  ,  (I)  ]  of  the  moving  observer  if  il>  and  d>m  are  related  by 
equation  (5).  This  leads  to  the  following  relation  : 


J  $m(wm)-du>m  *  /  $>st(u)-du>  (VuG  [0,®))  (A-l) 


m'  nr  m 

w  =0  u>=0 


Differentiating  equation  (5)  gives 


dw 


dw) 


m 


2U 

1  + - u> 

g 


Substituting  equation  (A-2)  in  (A-l)  gives  : 


(A-  2 ) 


d>  <*>„ 

pfTl  pm  s  t 

!  (  uj  )  ■  dw>  - 

J  m N  nr  ro  J 


nr  nr  m 

cj=0  u>  —0 

m  m 


2U 

1  + - u> 

g 


dw 


(A-  3) 


Differentiating  both  sites  of  equation  (A-3)  with  respect  to  <I>m 
gives  : 


$>_  (u>  ) 
m  '  nr 


2U 

1  + - w 


(A-4) 


g 
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MAXIMUM  VALUE  OF  THE  PSD 

Determine  the  frequency  at  which  the  power  spectral  density  function  has 
its  maximum  value. 

The  PSD  function  is  : 


•3.36 


k2 


0.84 


$  (u>)  = 
w 


'3rd 


kl 


(B-l) 


This  function  has  its  maximum  if  the  derivative  of  equation  (B-l)  equals 
zero.  The  derivative  is  : 


•  5  ■  kl  u> 
- e 


k2 
~ 4 


kl 


k2 

4 


-4  •  k2 


w 


(B-2) 


or 


k2 

4 

w 


-  e 


-  5 • kl • w  -  4 ■ kl • k2 

To 


(B-  3) 


From  equation  (B-3)  follows  if  $  equals  zero  : 


-  5 • kl • w  -  4 ■ kl • k2  -  0 


(B-4) 


So  the  frequency  at  which  the  PSD  function  has  its  maximum  value  is 


^  -  1.28  /  7h3  d 

vmax  ->ra 


(B-  5) 
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PSD  FROM  N  TIME  SAMPLES 

The  calculation  of  the  power  spectral  density  (PSD)  from  a  time 
dependent  variable  u  of  which  N  equidistant  samples  are  given  is 

presented . 

The  fourier  transform  of  u(t)  is  F(w)  : 

oc 

F(w)  -  f  u(t)-e->t  dt  (C-l) 

-  OC 

the  inverse  transform  is 

00 

u  ( t )  =  l/2n-f  F(u)  •  e  j dw  (C-2) 

-  oo 

Because  u(t)  is  only  sampled  in  the  interval  [-T/2  ,  T/2)  equation  (C-l) 
is  approximated  by: 

T/2 

Y(u)  -  f  u(t)  e'jwt  dt  (C- 3) 

-T/2 

T  is  the  observation  time. 

If  u(t)  is  sampled  every  At  seconds,  the  sampled  variable  u(t)  equals  : 
N/2-1 

u(t)  -  Z  u(t)-6(t  -  n-At)-At  (C-A) 

n--N/2 


with  At-T/N. 

The  fourier  transform  of  u(t)  is  : 
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T/2  T/2  N/2-1 

FQ(w)  -  /  u(t)-e'jwt-dt  -  /  2u(t)*(t 

-T/2  -T/2  n— -N/2 


n  •  At)  •  At  ■  e' •  dt  - 


N/2-1 

-  At-S  u(n-At)e'ja,nAt  (C-5) 

n— N/2 

Bv  renumbering  the  samples  with  un  -  u( (n-N/2-1) • At) ,  n-1 . N,  such 

that  the  first  sample  is  u^ ,  equation  (C-5)  becomes  : 

F.  (w)  -  At-ejw(NAt/2) -2  u  e-j“>(n-l)At  (C-6) 

n-1 


with  u>m  =  2m(m-l)/T  -  2?r(m-l)/(NAt)  and  m  -  1 . N  follows  : 


N 


At.ej*(ni-1)  .2  u  ,e-j27r(m-l)  (n-l)/N 


fm  =  Fu<<V  -  - 

n=l 


(C-7) 


=  -(-l)m-At  S  u  ,e-j2^(ni-l)(n-l)/N 


n-1 


f  is  the  discrete  Fourier  transform  component  of  frequency  of  the 
series  : 

(  un  ;  n-1, . N  } 

The  component  f  can  be  written  as  : 


-(-l)m(a  +  j  bm)_,At 


(C-8) 


Using  Euler's  formula  (e^x  -  cos(x)  +  jsin(x)),  the  componenta  am  and 

b  become  : 
m 
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N 

am  =  S  un-cos[ (m-1) (n-l)2x/N] 
n-1 


N 

bm  -  -  E  uR- sin[  (in- 1)  (n-l)2w/N] 
n-1 


(C-9) 


As  can  readily  be  seen,  the  information  about  the  Fourier  transform  is 

contained  in  the  first  N/2  values  of  am,  bm  :  (am,  bm  ;  m-l , . N/2 ) , 

because  the  spectrum  of  a  band  limited  signal,  sampled  every  At  seconds 
vanishes  for  u  values  higher  than  »r/At.  From  a  transform  point  of  view 
it  is  evident  that  if  the  time  signal  has  N  real  samples,  so  has  the 
transformed  time  signal.  This  is  because  of  the  invertibility  property 
of  the  Fourier  transform. 

The  amplitude  spectrum  is  found  from  : 

I  I  -  J  <-„2  +  *>/>  <c-10> 

With  a  fast  Fourier  transform  routine  the  first  N/2  coefficients  can  be 
calculated.  To  determine  the  PSD,  the  autocorrelation  function  must 
first  be  calculated.  The  sample  autocorrelation  function  Rjj(n)  of  u, 
based  on  N  samples,  is  determined  with  : 

1  N-n 

RN(n)  =  -  •  E  uj+n-uj  (C’n) 

N-n  j-1  J 

with  n  -  0 . nmax  and  RN(-n)  -  RN(n). 

The  PSD  function  is  the  Fourier  transform  of  the  autocorrelation 
function.  When  the  autocorrelation  based  on  a  finite  number  of  samples 
is  used  in  the  Fourier  transform,  an  approximation  S^(w)  of  the  PSD  is 
obtained  (see  ref.  4)  : 
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2 ' nraax+^ 

SN(-)-S  <C-12) 

n-1 

With  a  fast  Fourier  transform  routine  the  approximated  PSD  function  can 
be  determined  using  equations  (C-10),  (C-ll)  and  (C-12) . 
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